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Abstract 

Programmed cell death - apoptosis is one of the most studied biological 
phenomenon of recent years. Apoptotic regulatory network contains several 
significant control points, including probably the most important one - Bcl-2 
apoptotic switch. There are two proposed hypotheses regarding its internal 
working - the indirect activation and direct activation models. Since these 
hypotheses form extreme poles of full continuum of intermediate models, we 
have constructed more general model with these two models as extreme cases. 

By studying relationship between model parameters and steady-state re- 
sponse ultrasensitivity we have found optimal interaction pattern which re- 
produces behavior of Bcl-2 apoptotic switch. Our results show, that stimulus- 
response ultrasensitivity is negatively related to spontaneous activation of 
Bcl-2 effectors - subgroup of Bcl-2 proteins. We found that ultrasensitivity 
requires effector's activation, mediated by another STibgroup of Bcl-2 proteins 
- activators. We have shown that the auto-activation of effectors forms ul- 
trasensitivity enhancing feedback loop, only if mediated by monomers, but 
not by oligomers. Robustness analysis revealed that interaction pattern pro- 
posed by direct activation hypothesis is able to conserve stimulus-response 
dependence and preserve ultrasensitivity despite large changes of its inter- 
nal parameters. This ability is strongly reduced as for the intermediate to 
indirect side of the models. 

Computer simulation of the more general model presented here suggest, 
that stimulus-response ultrasensitivity is an emergent property of the direct 
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activation model, that cannot originate within model of indirect activation. 
Introduction of indirect-model-specific interactions does not provide better 
explanation of Bcl-2 functioning compared to direct model. 
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1. Introduction 

Apoptosis is the most ubiquitous type of programmed cell death, that has 
been observed and studied for more than one century [I] . It is a physiological 
phenomenon common for all multicellular organisms, necessary to tissue and 
body genesis and homeostasis [21 E]- Defects of its regulation may cause 
numerous diseases, including cancer, autoimmunity and neurodegenerative 
disorders Oil]. Recent research suggests, that apoptosis is distinct process 
known not only by its characteristic morphology and genomic destruction, 
but also as a process of progressive cellular disassembly with remarkably 
high level of complexity and flexible, multifaceted regulation. Apoptosis 
may be initiated by triggering events from within the cell or from outside 
the cell [1] . Apoptosis signaling then proceeds through multiple independent 
pathways 

Most of the apoptotic signaling pathways converge to mitochondria and 
may cause rapid increase in mitochondrial membrane permeabilization [51 E] ■ 
Mitochondrial outer membrane permeabilization (MOMP) leads to release of 
several pro-apoptotic factors, including cytochrome c and Smac/DIABLO [71 
|8] to cytosol. Cytosolic cytochrome c then associates with pro-apoptotic fac- 
tors, such as adaptor protein Apaf-1, to form Apoptosome complex [U [9]. 
At the same time, Smac/DIABLO inhibits lAPs - the inhibitors of apop- 
tosis [HI E] • Thus, mitochondria play pivotal role in apoptosis signaling by 
integrating and sensing incoming apoptotic signals and eventually responding 
by permeabilization of their outer membrane. 

Growing evidence demonstrate that MOMP is regulated by proteins of 
Bcl-2 family [TOl |TT] , which includes both pro-apoptotic and anti-apoptotic 
members. Proteins of the bcl-2 family are categorized by their relation 
to the apoptotic process and according to their Bcl-2 homology (BH) do- 
mains in their a-helical regions [T2l [T3t [H] . Anti-apoptotic members such as 
Bcl-2 itself, Bcl-xL and Bcl-w exhjibit four BH domains (BHl-4) [T^, pre- 
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vent MOMP and protect cells from a wide range of cytotoxic impacts [T3] . 
Pro-apoptotic members lack the BH4 domain and can be further divided to 
BH3-only proteins, such as Bid, Bik, Bim, Bad, Noxa and PUMA and mul- 
tidomain proteins, such as Bax, Bak and Bok [131 E]- BH3-only proteins 
can be activated by multiple pro-apoptotic signals and cytotoxic conditions, 
including cytokine deprivation or DNA damage [12]. Multidomain Bcl-2 
proteins, also termed 'effectors', once activated, permeabilize mitochondrial 
outer membrane (MOM) by allowing formation of oligomerized pores (MAC 
- mitochondrial apoptosis-induced channel) tl6l[T7|. The interplay between 
the three groups of Bcl-2 family determines MOMP commitment and forms 
so-called Bcl-2 apoptotic switch [ISl UHl [IB US] • 

Although there is general agreement that Bcl-2 proteins regulate MOMP, 
details of the Bcl-2 apoptotic switch mechanism remain controversial. It 
is still not clear, how various members of the Bcl-2 family interact between 
each other, and which of these interactions are decisive for MOMP induction. 
In particular, it is not resolved, whether the BH3-only proteins are able to 
activate multidomain effectors such as Bax and Bak directly, or whether 
they act indirectly, through neutralization of Bcl-2 anti-apoptotic sentinels 
to initiate MOMP [T5l 1201 [19] . The mechanism, based on the neutralization 
of the anti-apoptotic proteins, is known as indirect model [121 [HI [22] . 
Another, based on direct activation of effectors Bax and Bak by BH3 only 
proteins corresponds to the direct model [TUl [211 [22] • Another closely related 
question, deals with the manner by which anti-apoptotic proteins inhibit 
activity of effectors [151 [20l [I9]. Do the anti-apoptotic proteins bind non- 
activated Bax/Bak to prevent their activation, or do they bind and neutralize 
activated Bax/Bak? 

In the indirect model (see Fig. 1, Indirect model), effectors are activated 
spontaneously by conformational changes without presence of external fac- 
tors [231 [19]. Ill indirect model, the anti-apoptotic Bcl-2 proteins are bound 
to effectors blocking their spontaneous activation [23]. Signals for MOMP 
commitment activate BH3-only proteins, which can bind anti-apoptotic Bcl- 
2 proteins. Active BH3-only proteins can displace effectors from their anti- 
apoptotic relatives, thus promoting MOMP [251 [211 [231 [21]- Nevertheless, 
spontaneous activation of Bax and Bak, suggested in the indirect model, still 
lack convincing experimental confirmation j,26, ,19j. 

In the direct model (see Fig. 1, Direct model), effectors are activated 
solely by interactions with some members of the BH3-only proteins [23 [2H1 
[211 [301131]. Such members (e.g. Bid and Bim) are therefore called activators. 
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Other members of BH3-only proteins (such as Bad and Bik), are called 'en- 
ablers', 'de- repressors' or 'sensitizers', because they can bind anti-apoptotic 
Bcl-2 proteins, but do not activate Bax or Bak directly. Thus, enablers can 
help to promote MOMP by displacing the activators from anti-apoptotic 
Bcl-2 proteins [21] ■ 

Cellular switches such as Bcl-2 apoptotic switch are molecular mecha- 
nisms converting continuous incoming signals to two mutually distinct out- 
puts. Their role is to ensure unambiguous transitions between two differ- 
ent cellular states |32|. One of the necessary requirements to generate this 
behavior is ultrasensitive reaction mechanism IMf 13^ . Systems with 
ultrasensitive responses are define as systems which are more sensitive to 
stimulus changes than hyperbolic systems (systems with response given by 
the Michaelis-Menten equation) [351 ESI [32] . 

Comparison of the indirect and the direct model of Bcl-2 apoptotic switch, 
based on measuring of ultrasensitivity has already been done in work of Chen 
et al. [20] ■ Chen et al. concludes that the direct model is more plausible 
compared to its indirect alternative, since it leads to greater robustness, 
especially with respect to ultrasensitivity preservation. 

In the present work, we take more general approach. Out approach origi- 
nates from a different point of view on the Bcl-2 apoptotic switch controversy. 
Since interactions contained in the model of indirect activation may take 
their place in addition to those contained in the direct activation model [22] , 
hypotheses of indirect and direct activation are not entirely contradictory. 
Since models of indirect and direct activation can be considered as special 
cases of the more general continuum of models, we came out with the hybrid 
model that merges the interactions from both indirect and direct models. 
We used the hybrid model to browse space between hypotheses of indirect 
and direct activation, looking for parameter settings and interaction patterns 
under which model would preserve expected behavior. 

The idea of combining indirect and direct activation into hybrid model 
is not entirely new, its has been previously discussed [2S1 122] • However, it 
has never been examined whether such hybrid model could provide a better 
explanation of Bcl-2 apoptotic switch functioning compared to the models 
of indirect/direct activation, nor whether any new system properties may 
emerge from combining of those models. We have tried to answer this ques- 
tions here, as well as to comprehensively investigate importance of individual 
interactions between particular Bcl-2 family for the proper switch function- 
ing. 
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1.1. Mathematical model and its biological relevance 

We constructed hybrid model of biological switch formed by interactions 
of Bcl-2 family of proteins (scheme of this model is depicted in Fig. 2). 
Since the hybrid model combines the models of indirect activation and direct 
activation, by setting particular reaction rates to zero (effectively omitting 
corresponding interactions) we can reduce hybrid model to indirect, or di- 
rect one. To reduce the overall complexity of the generalized model, while 
preserve fundamental features of Bcl-2 family proteins, our models include 
only individual members representing whole class of Bcl-2 proteins with sim- 
ilar properties. Anti-apoptotic Bcl-2 family members, such as Bcl-2, Bcl-w, 
Bcl-xL and Mcl-l [2T1 [T3] are represented by Bcl-2. Moreover, pro-apoptotic 
Bcl-2 family members are grouped with respect to their composition of BH 
domains. The members containing multiple BH domains - effectors, for ex- 
ample Bax and Bak [151 El EI], are represented by Bax. BH3-only members 
of Bcl-2 family - Bid, Bim, Bad, Bik, Noxa, PUMA and others, are all as- 
sumed to be able to interact with effector proteins and thus induce their 
activation. Therefore they are represented by Act (Activators). 

Biological mechanism formed by Bcl-2 family of proteins processes multi- 
tude input signals into a single output signal - permeability of mitochondrial 
outer membrane. Activation of effectors leads to mitochondrial outer mem- 
brane permeabilization [71 EZl EH] • Among incoming signals are truncation 
of Bid by caspase-8 and activation of other BH3-only proteins by diverse cy- 
totoxic conditions, for example cytokine deprivation or DNA damage [T3l 123] . 
These input signals are represented in our model by single input stimulus - 
E, reaction rate of Act conversion to Act-a (reaction - 1). 

E 

1. Act — )■ Act-a 

Inhibition of activation of effector proteins Bax by anti-apoptotic Bcl-2 pro- 
teins is essential element of both models of Bcl-2 apoptotic switch. While the 
indirect model implies binding of non-activated Bax by anti-apoptotic Bcl-2 
proteins (reaction - 2) [391 El], the direct activation hypothesis assumes that 
Bax activity is inhibited solely by reversible binding and neutralization of 
activated Bax by anti-apoptotic Bcl-2 proteins (reaction - 3). 

2. Bcl-2 + Bax ^ Bcl-2~Bax 

3. Bcl-2 + Bax-a ^ Bcl-2~Bax-a 

Activated Act-a also reversibly bind anti-apoptotic Bcl-2, leading to their 
mutual neutralization (reaction - 4) [391 131] . 
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4. Bcl-2 + Act-a ^ Bcl-2~Act-a 



Moreover, Act-a is able to displace effectors Bax (indirect model) (reaction 
- 5), Bax-a (direct model) (reaction - 6) or both (hybrid model) from their 
complexes with Bcl-2 [211 123]. In the hybrid model, both Bax and Bax-a 
can compete to form the complex with Bcl-2 (reaction - 7). 

5. Bcl-2~Act-a + Bax ^ Bcl-2~Bax + Act-a 

6. Bcl-2~Act-a + Bax-a ^ Bcl-2~Bax-a + Act-a 

7. Bcl-2~Bax-a + Bax ^ Bcl-2~Bax + Bax-a 

Where the indirect model includes spontaneous activation of effector (reac- 
tion - 8), the direct model assumes that effectors are activated by Act-a 
(reaction - 9) [2HI EHl [301 El]- Both of these reactions are incorporated in 
the hybrid model. Spontaneous neutralization of Bax-a (reaction - 10) cor- 
responds to unmediated conformational changes that suppres its activity. 

8. Bax — 7- Bax-a 

9. Act-a + Bax — > Act-a + Bax-a 

10. Bax-a — )■ Bax 

The activity of effectors results in formation of mitochondrial apoptosis- 
induced channels. In accordance with conclusions of Martinez et al. [ID] 
we have modeled assembly of these channels as incremental growth of homo- 
oligomers constituted from monomeric effector units (reactions 11 and 12). 
Since MACs typically contain around 9 monomers, for simplicity we delimit 
oligomers growth up to 20 monomer units. We have verified that further 
growth of oligomer size does not qualitatively influence results proposed in 
this work (data not shown). 

11. Bax-a + Bax-a ^ MAC2 

12. Bax-a + MAC, ^ MAQ+i 

Additionally, all of the present compounds are continuously degraded. Bcl-2, 
Bax and Act are at the same time produced to balance their degradation 
flows. 

13. (All) ^ 

14. ^ Bcl-2 

15. -> Bax 

16. Act 
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1.2. Initial concentration and reaction rates 

Initial concentrations of particular species have been estimated in accor- 
dance with experimentally obtained data published in literature [TUl HU |121 
|20] . Concentrations of anti-apoptotic Bcl-2 and effector proteins (represented 
by Bcl2 and Bax, respectively) are considered to be in range of hundreds of 
nanomols, concentration of effectors was reported to be two times higher 
than that of Bcl-2 [TOl [20] . Concentration of BH3-only proteins (represented 
by Act) is considered to be much lower, in range of 1-20 nanomols [HI |32] 
(Table 1). 

The rates of binding of activators by anti-apoptotic Bcl-2 proteins has 
been estimated according to corresponding dissociation constants Kd pub- 
lished in literature |l3l HI], assuming reverse reaction rate 10~^s~^. The 
rest of reaction rates (binding of activated and non-activated Bax by anti- 
apoptotic Bcl-2 proteins, direct activation of Bax and its inactivation) have 
been chosen in accordance with previously published models |20l US, [19] . All 
the species are continuously degraded with degradation rate corresponding to 
half-life ti/2 = 180mm. Productions of inactive Bax, Act and Bcl-2 proteins 
are modeled by zero-order reactions parametrized to balance degradation 
under initial conditions (Tables 2 and 3). 

We used values of input stimulus E from lO^'* to lO^^mm^^. These 
values correspond to caspase-8 catalytic activity in Bid truncation, referred 
as ~ 10^M~^-s~^ |l6], assuming the number of active caspase-8 in range from 
10" - 10^ molecules per cell volume |47j . 

1.3. Implementation details 

All models were expressed in the CMDL (chemical model definition lan- 
guage) as well as SBML (the systems biology markup language) format 
(SBML format is available in supplementary). All the simulations were done 
by using CMDL/SBML ODE solver - ODEtoJava-dopr54-adaptive within 
Dizzy - Chemical kinetics simulations software [48j . 

2. Results 

In this work we focused on ultrasensitivity of Bcl-2 apoptotic switch and 
its dependence on particular interactions between the members of Bcl-2 fam- 
ily. To quantify the sensitivity of particular models, we used approach pro- 
posed in the work of Legewie et al. [36] (brief summary can be found in 
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Appendix). As the reference response, used to calculate the relative amplifi- 
cation coefficient, we used Michaelis-Menten (hyperbolic) equation. Models 
with parameter sets, exhibiting relative amplification coefficient tir higher 
than unity were classified as ultrasensitive [36] . 

First of all, we compared steady-state stimulus-response dependence (Fig. 
3) of hybrid, indirect and direct models (using reference setup, parameters 
from Table 3). The input stimulus is represented by the value of the param- 
eter E - rate of Act to Act-a conversion. 

In most of the previously published modeling works response have been 
quantified by the relative activity of effectors. In contrast, in this work 
we quantified response as the sum of amounts of individual oligomers that 
contains more that 6 monomer units, weighted by their size (Eq. 1). Such 
calculation is reflecting the fact that MAC pore complex needs to contain 
a minimum 6 Bax/Bak monomer units to allow transport of cytochrome 
c (the major pro-apoptotic factor released from mitochondria), also taking 
into account that the larger pores contribute to permeabilization more then 
smaller ones [40j. 

20 

Response = ■ MAQ (1) 

i=6 

For all three models we evaluated the dependence of the response coefficient 
on the stimulus, and the activated fraction [IHl EO]- Then we calculated 
corresponding relative amplification coefficients nn. 

While the indirect and hybrid models were proved to be strongly sub- 
sensitive {nn ~ 10~^), the direct model resembled ultrasensitive behavior 
{riR = 2.0). 

2.1. Variations of pivotal parameters and their effects on sensitivity 

In order to investigate the influence of individual reaction parameters on 
sensitivity of the Bcl-2 apoptotic switch models, we plotted relative ampli- 
fication coefficient of hybrid model as a function of variation of individual 
reaction parameters (ki, kil, ki2, kin, kc, ks, ko). The same parameter vari- 
ation was then performed after the reduction of the hybrid model into its 
limiting direct /indirect cases. 

Our results demonstrate, how dramatically can variation of single reac- 
tion parameter effect the sensitivity of the modeled system (Fig. 4). The 
variation of the parameter ks (spontaneous activation of Bax) change the rel- 
ative amplification coefficient of hybrid model by two orders of magnitude. 
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As ks is lOOx decreased hybrid model reaches the ultrasensitivity. Less ap- 
parent is the influence of parameters kc (activation of Bax by Act-a) and 
kin (inactivation of Bax-a). The dependence of sensitivity to parameter's 
variation is monotonic for most of the parameters. One notable exception 
is non-monotonic dependence on reaction parameter ki (binding of Bax by 
Bcl-2). We found local maxima of sensitivity around ki — 12.0 • 10~^ and 
ki = 7.0 ■ 10~^, for the hybrid and indirect model respectively. Similarly, the 
dependence of relative ampliflcation coefficient of the direct model on varia- 
tion of parameters contains several local extreme values. By comparing our 
models, we can see how different is their capability to preserve ultrasensitive 
response. While direct model keeps ultrasensitive response under wide range 
of variations, indirect model yields response only with relative amplification 
coefficient far below unity. 

Although single point variation of reaction parameters may reveal certain 
interdependence of system properties on the value of chosen parameter, it 
provides only very limited view on relationships between these properties 
and system parameters. 

To get a deeper understanding of the impact of the individual Bcl-2 pro- 
teins interactions and their interplay on the switching behavior of the mod- 
eled system, we performed an additional analysis. We have generated 2000 
sets of randomly changed values of reaction rates ki, kil, kc and ks related 
to controversial interactions. At the same time the rest of the parameters 
remained unchanged at their reference values. Then we calculated relative 
amplification coefficient of the hybrid model corresponding to each set of 
parameters. To obtain picture about how the particular parameters affect 
sensitivity of hybrid model, we separated all the parameter sets to intervals, 
according to the variation of individual parameters. Then, for each interval 
we calculated the mean of relative amplification coefficients - ur given by the 
parameter sets within the interval. 

As can be seen from Fig. 5, the mean relative amplification coefficient of 
parameter sets, in which the value of ks is decreased (compared to its refer- 
ence value) , is much higher than the mean relative amplification coefficient of 
parameter sets in which the ks is increased. This means, that the decreased 
value of ks yield higher relative amplification coefficient more often than the 
increased ones. In case of the parameter kc, the highest mean relative am- 
plification coefficient (almost ur = 0.8) was obtained for parameter sets in 
which kc was varied from 17.8 to 31.6 x its reference value. There is no 
such apparent dependency observed in case of parameters ki and kil. The 
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mean relative amplification coefficient is approximately the same, regardless 
of variation of these parameters (data not shown). 

Our results indicate that the presence of spontaneous activation of ef- 
fectors (with corresponding parameter ks) adversely influence sensitivity of 
hybrid model. On the contrary, activation of effectors by activators (param- 
eter kc) seems to be beneficial for ultrasensitive, switching behavior of the 
model. These findings are unfavorable for the hypothesis of indirect activa- 
tion, that rely on the spontaneous activation of effectors without activation 
by activators. 

Taking into account the findings above, we set parameters ks and kc of 
hybrid model to 10~^x and 20.0 x their reference values, respectively. In next 
step we investigated how the interplay between two alternatives of effectors 
inhibition may affect sensitivity of the hybrid model. We have simultane- 
ously varied values of parameters ki and kil (corresponding to inhibition of 
effectors before activation, and after activation, respectively) and measured 
the relative amplification coefficient. 

Results shows (Fig. 6) that sensitivity of the hybrid model is only weakly 
affected by changing the values of ki and kil. Relative amplification coef- 
ficient of hybrid model is changed only marginally, even as we changed the 
values of ki and kil by several orders of magnitude. It clearly appears that 
the highest relative amplification coefficient is reached when the value of the 
parameter ki is reduced and of the parameter kil is increased, indicating the 
importance of neutralization of active effectors by anti-apoptotic proteins. 
The lowest relative amplification coefficient was found as we reduced values 
of both parameters. 

2.2. Addition of the auto- activation of Bax 

Works of Ruffolo et al. [50] and Tan et al. [51] provided experimental 
evidence that activated effectors Bak, as well as Bax can directly interact 
with, and positively influence activity of non-activated Bak and Bax, respec- 
tively. This interaction, called Bax/Bak auto-activation can enhance activity 
of effectors, initially activated by spontaneous activation (as is proposed by 
hypothesis of indirect activation) or by BH3-only activators (hypothesis of 
direct activation). 

Such auto-activation forms positive feedback loop, that can dramatically 
change the stimulus-response steady-state properties of the model. This 
could make our search for ultrasensitive setup even more difficult, therefore 
we have not involved this interaction into initial stages of our model. On 
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the other hand, addition of such positive feedback to ultrasensitive reaction 
mechanism may lead to increase of the ultrasensitivity. 

It is still not determined whether auto-activation of effectors is mediated 
by monomeric, or by oligomerized effectors |26]. In most of the previous 
modeling works, auto-activation of effectors was modeled as arrangement of 
the inactive effectors into oligomerized pores. In contrast to these works, 
we have tried to investigate, which of the two alternatives of auto-activation 
has the most beneficial effect on the switching properties of the Bcl-2 apop- 
totic switch. We have modeled the auto-activation mediated by monomeric 
effectors as the following reaction (reaction autol): 

autol: Bax + Bax-a Bax-a + Bax-a 

Auto-activation mediated through oligomerization was modeled as an ar- 
rangement of the inactive Bax into oligomerized pore (reaction auto2), sim- 
ilarly to previous models of Bcl-2 apoptotic switch [191 IS2] ■ 

auto2: Bax + MAC^ MAC^+i 

We have added auto-activation of Bax into the hybrid model and mea- 
sured relative amplification coefficient for different values of corresponding 
reaction rate - ka. We tested first the impact of the auto-activation of Bax 
mediated by monomeric units of active Bax-a. We have found that for narrow 
area of ka values rapid rise of the relative amplification coefficient is followed 
by dramatic decline as value the of ka exceeds the threshold ~ 2.0 ■ 10~^ (see 
Fig. 7). Increasing the value of the parameter ka amplifies the strength of the 
positive feedback. As we take a closer look on the stimulus-response curves 
(Fig. 8), we can see that the positive feedback narrows the intermediate part 
of the sigmoidal stimulus response curve, making it steeper. 

There are numerous works discussing how the addition of positive feed- 
back into an ultrasensitive reaction system may influence its steady-state 
properties [53l [Ml ESj. Such addition is often associated with emergence of 
the bistable behavior assuming the feedback strength exceeds certain thresh- 
old. In contrast with these suggestions, we observe rapid fall of the relative 
amplification coefficient for overly strong feedback. The inflection point of 
the stimulus-response curve shifts to the left along the stimulus axis as we 
increase the ka value, and when the certain threshold value is reached, the 
inflection point leaves reasonable range of input stimuli. The auto-activation 
then causes high mitochondrial membrane permeability even at basal level 
of incoming stimuli. 
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We have replaced auto-activation mediated by monomers - reaction autol, 
by reaction auto2 and evaluated the relative amplification coefficient for dif- 
ferent values of corresponding reaction parameter - ka. Relative amplification 
coefficient was only slightly affected by variation of the ka value, indicating 
that auto-activation mediated by oligomerized effectors has only marginal 
effect on the sensitivity of the hybrid model. 

These results show, that although both forms of auto-activation may 
efficiently enhance the effectors activity, only the auto-activation mediated 
by monomers interaction can amplify the steady-state ultrasensitivity of the 
hybrid model. Although, it should be noted that sensitivity increase occurs 
only withing the narrow range of values of auto-activation reaction rate. 
Overly fast auto-activation of effectors causes rapid onset of permeability 
even if the input stimulus is weak, thus disabling switching behavior. 

2.3. Robustness analysis 

All essential biological mechanisms should exhibit considerable level of 
robustness against parameter changes |56] , since effects of environmental dif- 
ferences, polymorphism or mutations needs to be compensated [57]. Biolog- 
ical switches such as Bcl-2 apoptotic switch should preserve ultrasensitivity 
as their pivotal feature, in spite of variation of parameters. 

First, we have adjusted the parameter setup of the hybrid model in ac- 
cordance to all the previous results, to obtain desirable level of sensitivity 
(adjusted values of parameters are summarized in Table 4). To test the 
robustness of our models, we measured relative amplification coefficient for 
multiple sets of varied reaction rates (ki, kil, ki2, kc, ks, kin) and initial con- 
centrations (Act, Bax, Bcl-2). Similar to the work of Barkai and Leibler [5H] 
and some later works [571 EO], we generated set of modified parameters, 
by multiplying each of reference parameters by 10'', where q is random real 
number taken from Gaussian distribution (mean = 0.0, variance = 1.0). Rel- 
ative amplification coefficient of each set was then plotted as a function of 
total parameter variation (T), which is defined as total order of magnitude 
of parameter variation [5HI ETl 120] 



where rip is number of varied parameters, Pi is varied parameter and Pi^ref is 
corresponding reference parameter. Robustness was quantified by frequency 
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of occurrence of ultrasensitive responses among 1000 random parameter sets. 
Similar approach was previously used in the work of Chen et al. [20|, where 
direct and indirect model were compared. It needs to be remarked that 
our model considerably differs to those proposed by Chen et al., therefore 
robustness analysis results are not directly comparable. 

We have found around 35% from 1000 hybrid model parameter sets which 
yield sigmoidal response curve with corresponding relative amplification co- 
efficient higher than 1.0 (see the scatter plot in Fig. 9). The values of the 
relative response coefficient mostly occurred below the 2.5, in around 90% 
of all sets. As we reduced the hybrid model to the indirect model, none 
of the thousand sets yield sigmoidal response curve with relative amplifi- 
cation higher than 1.0, strongly opposing the suggestion, that the indirect 
model gives rise to ultrasensitive response. On the other hand, as we reduced 
the hybrid model to the direct model, we have obtained sigmoidal response 
curve with relative amplification coefficient higher that one in around 60% 
of parameter sets. This shows that omitting the spontaneous activation of 
effectors, and their preventive inhibition makes response more model more 
sensitive. 

Frequency of occurrence of relative amplification coefficient higher then 
unity, is a measure of model's capability to preserve ultrasensitive stimulus- 
response dependence despite changing its internal parameters. Although 
ability to preserve ultrasensitive response is directly related to switching ro- 
bustness, there is something else we would expect from the robust ultrasen- 
sitive switch. In fact, besides preserving ultrasensitivity of stimulus- response 
dependence, robust switch should be able to conserve the stimulus-response 
dependence itself. With respect to subsequent processes regulated by ultra- 
sensitive switch, we consider, that proper functioning of such switch requires 
reaction mechanism that is able to conserve its stimulus-response character- 
istics despite small changes of its internal parameters. We hypothesized, that 
spontaneous activation of effectors, or the preventive inhibition of effectors 
may possibly strengthen this capability. To characterize the model's ability 
to preserve the stimulus-response dependence, we have measured the position 
of its inflection point. 

We took parameter sets, which yield relative amplification coefficient 
higher than one, and plotted the corresponding stimulus-response curves and 
stimulus- response slope dependencies (see Fig. 10). For each response de- 
pendence, we localized the inflection point at which response slope reaches 
its maximum. Then, we plotted statistical distribution of inflection point 
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coordinates (Fig. 11). As can be seen from Fig. 11, for the direct model, 

there arc two statistically significant peaks located around stimulus values 
1.8 ■ 10^'^ and 9.5 ■ 10^"^. This indicates, that the interaction pattern of di- 
rect model is efficiently conserving stimulus coordinate of the inflection point 
in close proximity of these values. Hereby, stimulus coordinate of inflection 
point changes only a little as we change internal parameters of the direct 
model. We do not observe such conservation of stimulus coordinate in the 
case of the hybrid model (Fig. 11 - Hybrid model), showing that stimulus 
coordinate of inflection point varies proportionally to variation of internal 
parameters. 

The distributions of the response coordinate of the inflection point of both 
models show, that the inflection point occurs typically around 0.1 to 0.2 % 
of the highest response reached. The response coordinate of the inflection 
point is strongly conserved in both hybrid/direct models. 

Taken together, we found that the direct model is, compared to the hy- 
brid one, more robust with respect to ultrasensitivity, as well as regarding 
the switching threshold behavior of inflection point of its stimulus-response 
dependence. 

3. Discussion 

3.1. Steady-state ultrasensitivity originates from Bax/Bak activation by BH3- 
only activators, not from spontaneous activation 

Our aim was to identify those interactions between Bcl-2 subgroups, 
which are responsible for emerging switching behavior. Therefore we ex- 
amined how the variations of the model's reaction rates can affect properties 
in the modeled system. Single parameter variations, as well as simultaneous 
variations of multiple parameters revealed interesting relations between pa- 
rameter setup and steady-state sensitivity of the hybrid model and reduced 
indirect/direct models. 

Based on these results we can conclude that spontaneous activation of ef- 
fectors Bax/Bak plays minor if any role in activation, allowing their further 
oligomer izat ion and mitochondrial permeabilization. Similarly, it seems that 
the neutralization of activated effectors contributes to switching behavior 
far more than preventive inhibition of inactive forms. Taken together, the 
interaction pattern corresponding to hypothesis of direct activation forms 
switching mechanisms that yields ultrasensitive behavior. Alternative inter- 
action pattern - the indirect model, does not lead to this behavior. 
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3.2. Monomers mediated auto- activation of Bax/Bak may strongly enhance 
ultras ensitivity of the Bcl-2 apoptotic switch 

As was reported in work of Tan et al. [51], activated Bax may activate 
non-activated Bax. Similar finding was reported for Bak in work of Ruffolo 
et al. [5^. The mechanism of this auto-activation is still not fully under- 
stood. It is discussed, whether auto-activation is mediated by monomers or 
by oligomerized Bax/Bak. Despite the lack of knowledge, auto-activation of 
effectors is sometimes referred as responsible for switching properties of Bcl-2 
apoptotic switch. The relevance of the auto-activation with respect to the 
switching function of Bcl-2 regulatory mechanism has never been examined. 

Such autocatalytic activations of effectors form positive feedback loop, 
that can strengthen ultrasensitivity of the reaction mechanism or even lead 
to bistable behavior [M]. We have modeled the influence of auto-activation 
on the switching properties of the studied system, by measuring sensitivity 
after addition of auto-activation and variation its reaction rate. 

We have added autocatalytic activation of effectors into adjusted model 
and varied its reaction rate. In case of monomers mediated auto-activation, 
within certain values of reaction rate we observed rapid increase of sensi- 
tivity, followed by its sudden fall. This is due to the withdrawal of the 
inflection point of sigmoid stimulus-response curve for the reasonable range 
of input stimuli. This means that overly strong auto-activation may easily 
cause massive effectors activation and high MOM permeability even by basal 
stimuli. Auto-activation mediated by oligomers has no significant effect on 
the ultrasensitivity of the modeled system. 

Although hybrid model of the Bcl-2 apoptotic switch exhibits sufficient 
level of ultrasensitivity even without presence of any auto-activation feedback 
loop, addition of monomers-mediated auto-activation may strongly amplify 
ultrasensitivity. Ultrasensitivity amplification is attributed to monomers me- 
diated auto-activation only. Auto-activation by oligomerized effectors can 
only enhance activity of effectors but does not contribute to ultrasensitive 
behavior. 

3.3. The direct model is more robust, compared to the hybrid or indirect ones 
Upon addition of the auto-activation into the hybrid model, we have ad- 
justed reaction rates of this model to obtain strongly ultrasensitive response. 
Then we analyzed its robustness with respect to its ability to preserve ultra- 
sensitivity against variations of its reaction parameters and initial conditions. 
We measured frequency of occurrence of the ultrasensitive response among 
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thousand of randomly generated parameter sets. There were around 35 % of 
the parameter sets of the hybrid model which yielded sigmoidal stimulus re- 
sponse with relative amplification coefficient higher than one. We have then 
reduced the hybrid model to the indirect/direct ones and compared robust- 
ness of the hybrid model with robustness of indirect and direct models. We 
have found that indirect model is unable to provide ultrasensitive behavior, 
since none of the thousand parameter sets have yield sigmoidal stimulus re- 
sponse dependence with relative amplification coefficient higher than one. In 
opposite, the direct model compared to the hybrid model was found much 
more robust (70 % of the parameter sets), indicating that the model's abihty 
to preserve ultrasensitivity is negatively affected by spontaneous activation 
of effectors and their preventive inhibition. 

We have hypothesized that these interactions may possibly increase model's 
ability to conserve stimulus-response dependence. In order to examine this 
assumption, we have compared stimulus response curves of ultrasensitive re- 
sponses obtained from previous analysis. We have plotted the distribution 
of the threshold stimulus (the stimulus corresponding to inflection point of 
the stimulus response sigmoid) and the corresponding threshold response. 
In case of the direct model, we have observed two significant peaks of that 
distribution. This interesting finding indicates that the direct model exhibits 
the capability to conserve threshold stimulus within area around those peaks. 
We found no such effect in case of the hybrid model. This rules out our as- 
sumption about robustness supporting role of the spontaneous activation and 
preventive inhibition of effectors. 

3.4- Indirect model interaction may provide some other system properties, 
but not steady-state ultrasensitivity 

Taken together, all the results from this work point to following conclu- 
sions: interaction pattern as described by the hypothesis of indirect activa- 
tion is very unlikely to give arise to the switching behavior as is expected 
from Bcl-2 apoptotic switch. On the other hand, its alternative - the direct 
model hypothesis gives much more plausible explanation of functioning of 
Bcl-2 apoptotic switch. By merging these interaction patterns into one, we 
obtained hybrid model. Hybrid model, compared to direct one, is less robust 
with respect to ultrasensitivity preserving, and seems to be less efficient in 
providing switching function. 

Despite these conclusions related to steady-state ultrasensitivity, we can- 
not completely rule out the indirect model interactions. Since the require- 



16 



ments for ultrasensitivty suggest that these interactions must proceed at very 
low rates, but these interactions still may exist and provide other important 
system properties. Low rate spontaneous activation of effectors can make 
response to changed stimuli more abrupt, while still preserving desirable 
steady-state ultrasensitivity. At certain affinity, continuous binding and in- 
hibition of inactive effectors by anti-apoptotic proteins may possibly help to 
attenuate the noise of the incoming signal. However these ideas need futher 
examination. Therefore, the study of transient system properties of the Bcl-2 
apoptotic switch is aim of our further modeling work. 
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5. Appendix 



In metabolic control analysis, the response coefficient is defined as steady 
state fractional change of the system response AX/X divided by fractional 
change of stimulus AS/ S in limit as AS tends to zero 



R 



x 



AX/X _ dlnX 
A^^^o AS/S ~ dlnS 



lim 



(A.l) 



The global sensitivity is measured through relative amplification coefficient 
using definition of Legewie et al. |36j as: 



R^df 



riR 



fE 



(A.2) 
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Where / is activated fraction defined as: 

/ = / (A.3) 

fu and fi are margins of the activated fraction range and R^^ is the response 
coefficient of the reference response X^, which can be any monotonically 
increasing or decreasing function jSH] • 
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Figure 1: Layout of the interplay between particular members of the Bcl-2 family as 
described by the indirect activation model (top) and direct activation model (bottom). 
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Figure 2: Simplified scheme of the hybrid model as proposed in this work (production and 
degradation flows, as well as reverse reaction are not depicted). 
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Figure 3: Plot of the stimulus-response dependence (response is scaled to one) of hybrid 
model, as well as indirect and direct models, using reference parameter setup (see Table 
3). 




Figure 4: Dependence of relative amplification coefficient of indirect hybrid (left) and 
reduced - indirect (center), direct (right) models, on variation of individual parameters. 
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Figure 5: Dependence of the mean relative amplification coefficient on variation of 
eter as used in randomly chosen parameter sets. 
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Figure 6: Dependence of relative amplification coefficient on simultaneous variation of 
reaction parameters. Parameters ki/kil are reaction rates of binding and inhibition of 
inactive/active effectors by their anti-apoptotic relatives. 




Figure 7: Dependence of relative amplification coefficient on the variation of reaction rate 
of Bax auto-activation. Black line corresponds to auto-activation mediated by monomers 
(as described by reaction autol), gray line corresponds to auto-activation mediated by 
oligomerized Bax. 
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Figure 8: Changing shape of stimulus-response curves upon addition of monomers me- 
diated Bax auto-activation (autol) with increasing reaction rate (ka values are displayed 
beside each curve). 
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Figure 9: Scatter plots of relative amplification coefficient with dependence on the value of 
total parameter variation. Dots localized over gray dashed line corresponds to randomly 
generated parameter sets that yield ultrasensitive response, with relative amplification 
coefficient higher than one. 
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Figure 10: For each randomly generated parameter set yielding ultrasensitive response 
(see Fig. 9, dots over gray dashed line), corresponding stimulus-response dependence 
was plotted (gray lines). Stimulus-response dependencies arc normalized, and merged to 
common inflection point. Similarly, corresponding stimulus-response slope dependencies 
are plotted. Black lines correspond to response produced under adjusted parameter setup 
(see Table 4). 
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Figure 11: For each stimulus-response dependence (Fig. 10, gray lines) we have measured 
its threshold stimulus and corresponding threshold response. Histograms shows distribu- 
tion of these values. 
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Species Initial concentration Notes & Ref. 





# • cell-^ 


nM 




Bcl-2 


60 000 


100 


[iniEo] 


Bax 


120 000 


200 


[I0l[20] 


Act 


6000 


10 


[111112] 


Bax-a 










Act-a 










Bcl-2~Bax 










Bcl-2~Bax-a 










Bcl-2~Act-a 










MACi 











Table 1: List of initial concentrations. Concentrations are listed as numbers of molecules 
per reference cell volume 1 pi, and in more common units - nanomols. 



No. 


Reaction 






1. 


E 

Act — > Act-a 


E 




2. 


Bcl-2 + Bax Bcl-2~Bax 


ki 


kmi 


3. 


Bcl-2 + Bax-a Bcl-2~Bax-a 


kil 


kmil 


4. 


Bcl-2 + Act-a Bcl-2~Act-a 


ki2 


kmi2 


5. 


Bcl-2~Act-a + Bax Bcl-2~Bax + Act-a 


ki 


ki2 


6. 


Bcl-2~Act-a + Bax-a Bcl-2~Bax-a + Act-a 


kil 


ki2 


7. 


Bcl-2~Bax-a + Bax Bcl-2~Bax + Bax-a 


ki 


kil 


8. 


Bax —7- Bax-a 


ks 




9. 


Act-a + Bax — t- Act-a + Bax-a 


kc 




10. 


Bax-a — )• Bax 


kin 




11. 


Bax-a + Bax-a MAC2 


ko 


kmo 


12. 


Bax-a + MACi ^ MACj+i 


ko 


kmo 


13. 


(All) ^ 


kd 




14. 


Bcl-2 


kpBcl-2 




15. 


— )• Bax 


kpBax 




16. 


Act 


kpAct 





Table 2: List of reactions and corresponding reaction parameters. Reversible reactions are 
listed with forward/reverse reaction parameters. 
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Par am. 



Value 



Unit 



Notes & Ref. 



ki 


1 


0-10-6 


(1.0-10^^) 


ce//-#- 


-1 


min~ 


-1 


(M-i 


.-1) 


(in direct) = 0, estimated 


kmi 




0.06 


(0.001) 






min" 


-1 


(s-i) 




= kmil 


kil 


1 


0-10-6 


(1.0-10^) 


cell-#' 


-1 


min" 


-1 


(M-i 




(in indirect) = 0, estimated 


kmil 




0.06 


(0.001) 






min^ 


-1 


(.-1) 




same as in [TH] 


ki2 


1 


0-10-6 


(1.0-10^^) 


cell-#- 


-1 


min~ 


-1 


(M-i 




Kd = lOOnM 031 m 


kmi 2 




0.06 


(0.001) 






min" 


-1 


(.-1) 




same as in [19j 


kc 


1 


0-10-6 


(1.0-10"^) 


cell-#~ 


-1 


min" 


-1 


(M-i 




(in indirect) = 0, 


ks 




0.06 


(0.001) 






min' 


-1 


(.-1) 




estimated 


kin 




0.06 


(0.001) 






min~ 


^1 


(.-1) 




estimated 


ko 


1 


0-10-6 


(1.0-10^^) 


cell-#~ 




min~ 


^1 


(M-i 




estimated 


kmo 




0.06 


(0.001) 






min~ 


^1 


(.-1) 




estimated 


kd 




0.0039 


(6.5-10-5) 






min^ 


^1 


(.-1) 




ti/2 = 180 min, estimated 


kpBcl-2 




234.0 


(6.4-10-12) 


#-cell- 


-1 


min^ 


-1 


{M-s' 




= Bcl-2j„irkd 


kpBax 




468.0 


(1.3-10-11) 


#-cell- 


-1 


min' 


-1 


{M-s- 




= Baxjnirkd 


kpAct 




23.4 


(6.4-10-13) 


#-cell~ 


-1 


min~ 


-1 


{M-s- 







Table 3: List of reference values of reaction parameters. By setting parameters kil and 
kc equal to zero, the hybrid model is reduced to the indirect one. Similarly, the direct 
model can be obtained by setting parameters ki and ks of the hybrid model equal to 
zero. Production rates kpBcl-2, kpBax and kpAct has been set to balance degradation of 
corresponding species under initial conditions. 



Par am. Value Unit 



ki 


1.0 


10" 


-8 


(1.0 


102) 


cell-#- 


-1 


min~ 


^1 


(M-i 


s' 




kil 


1.0 


10" 


-4 


(1.0 


106) 


cell-#- 


^1 


min' 


-1 


(M-i 


s' 




kc 


2.0 


10" 


-5 


(2.0 


105) 


cell-#' 


^1 


min' 


^1 


(M-i 


s' 




ka 


1.8 


10" 


-6 


(1.8 


10^) 


cell-#- 


-1 


min~ 


-1 


(M-i 


s' 




ks 


6.0 


10" 


-4 


(1.0 


10-5) 






min~ 


-1 


(.-1) 







Table 4: List of adjusted values of reaction parameters. In order to obtain desirable level of 
ultrasensitivity, listed parameters have been adjusted according to previous results. Rest 
of parameters remained unchanged. 
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